library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)

filekp=read.table("ASM_Bayes_Meta_Chisq_anno_num2.txt",head=T,sep="\t")
filekp[filekp=="."]=NA

sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel2=c("M30_M29","M26_M25","M35_M36","M28_M27")
#sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M42_M41","M44_M43","M48_M47","M50_M49","M52_M51");sel2=c("M28_M27","M30_M29","M26_M25","M36_M35")
sel3=c("M18_M17","M20_M19","M22_M21","M40_M39")

sel3=c("M18_M17","M20_M19","M22_M21","M40_M39")

sel1=unlist(strsplit(sel1,"_"))
sel2=unlist(strsplit(sel2,"_"))
sel3=unlist(strsplit(sel3,"_"))

sel=c(sel1[seq(1,length(sel1),2)])#con
sel=c(sel1[seq(2,length(sel1),2)])#case
sel=c(sel1[seq(2,length(sel1),2)],sel2)
sel=c(sel1[seq(1,length(sel1),2)],sel3)
sel=sel2
sel=sel3
#sts=c(rep(c(0,1),length(sel1)/2),rep(c(1,1),length(sel2)/2))
file=filekp[,c(paste(sel,"_reads1",sep=""),paste(sel,"_reads2",sep=""))]

result=matrix(,nrow(file),ncol=4)
for (i in 1:nrow(file)){
reads=as.numeric(as.matrix(file)[i,])
ref=reads[1:(ncol(file)/2)]
var=reads[(ncol(file)/2+1):ncol(file)]
tid=(1:(length(sel1)/2))-1
#tid=c((1:(length(sel1)/2))-1, rep((length(sel1)/2+1):((length(sel1)+length(sel2))/2),each=2)-1)
#tid=c((1:(length(sel1)/2))-1, rep((length(sel1)/2+1):((length(sel1)+length(sel3))/2),each=2)-1)
#tid=(1:(length(sel)))-1		#双发和健康单独为组的组号
sl=which(!is.na(ref))
if ((length(sl)/2)>=1){
#df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x1=sts[sl],x2=tid[sl])
df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x2=tid[sl])
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1 + f(x2, model = "iid")
#formula = ALT_COUNT ~ 1 + f(TISSUE_ID, model = "iid")
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}else{
#df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x1=sts[sl])
df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl])
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}
}
filekp$DC_con_pvalue=result[,1]
filekp$DC_con_mean=result[,2]
filekp$DC_con_down=result[,3]
filekp$DC_con_up=result[,4]

filekp$DC_case_pvalue=result[,1]
filekp$DC_case_mean=result[,2]
filekp$DC_case_down=result[,3]
filekp$DC_case_up=result[,4]

filekp$DC_case_CC_pvalue=a[,1]
filekp$DC_case_CC_mean=a[,2]
filekp$DC_case_CC_down=a[,3]
filekp$DC_case_CC_up=a[,4]

filekp$DC_con_HC_pvalue=a[,1]
filekp$DC_con_HC_mean=a[,2]
filekp$DC_con_HC_down=a[,3]
filekp$DC_con_HC_up=a[,4]

filekp$CC_pvalue=result[,1]
filekp$CC_mean=result[,2]
filekp$CC_down=result[,3]
filekp$CC_up=result[,4]

filekp$HC_pvalue=a[,1]
filekp$HC_mean=a[,2]
filekp$HC_down=a[,3]
filekp$HC_up=a[,4]

dropit=c("intergenic","ncRNA_intronic","ncRNA_exonic","ncRNA_splicing")
kk=which(!filekp$Func.refGene %in% dropit)
file1=filekp[kk,]
file1=file1[file1$DC_con_pvalue_fdr<=0.05 | file1$DC_case_pvalue_fdr<=0.05,]
file1=file1[file1$Bayes_SZBD_in_DC>=100,]
#########################以下为针对每个样本算ASE

filekp=read.table("ASM_Bayes_Meta_Chisq_anno_num2.txt",head=T,sep="\t")
filekp[filekp=="."]=NA
sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49");sel2=c("M28_M27","M30_M29","M26_M25","M36_M35")
sel3=c("M18_M17","M20_M19","M22_M21","M40_M39")
sel1=unlist(strsplit(sel1,"_"))
sel2=unlist(strsplit(sel2,"_"))
sel3=unlist(strsplit(sel3,"_"))
sel=c(sel1,sel2,sel3)
file=filekp[,c(paste(sel,"_reads1",sep=""),paste(sel,"_reads2",sep=""))]
result=data.frame(matrix(filekp$unitID,nrow(file),ncol=1))
names(result)=c("unitID")
for (i in 1:nrow(file)){
reads=as.numeric(as.matrix(file)[i,])
ref=reads[1:(ncol(file)/2)]
var=reads[(ncol(file)/2+1):ncol(file)]
sl=which(!is.na(ref))
for(k in 1:length(sl)){
df=data.frame(y=var[sl[k]],Ntrials=ref[sl[k]]+var[sl[k]])
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	name=names(file)
	sample_ase_p=paste(gsub("_reads1","",name[sl[k]]),"_ASE_Pvalue",sep="")
	sample_ase_mean=paste(gsub("_reads1","",name[sl[k]]),"_ASE_mean",sep="")
	result[i,sample_ase_p]=post_pred_p
	result[i,sample_ase_mean]=coef$mean
	}
}

###########################以下为针对健康和双发分别为组算ASE
for (i in 1:nrow(file)){
reads=as.numeric(as.matrix(file)[i,])
ref=reads[1:(ncol(file)/2)]
var=reads[(ncol(file)/2+1):ncol(file)]
#tid=(1:(length(sel1)/2))-1
#tid=c((1:(length(sel1)/2))-1, rep((length(sel1)/2+1):((length(sel1)+length(sel2))/2),each=2)-1)
#tid=c((1:(length(sel1)/2))-1, rep((length(sel1)/2+1):((length(sel1)+length(sel3))/2),each=2)-1)
tid=(1:(length(sel)))-1		#双发和健康单独为组的组号
sl=which(!is.na(ref))
if ((length(sl)/2)>=1){
#df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x1=sts[sl],x2=tid[sl])
df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x2=tid[sl])
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1 + f(x2, model = "iid")
#formula = ALT_COUNT ~ 1 + f(TISSUE_ID, model = "iid")
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}else if(length(sl)>=1){
#df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x1=sts[sl])
df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl])
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}
}

#######################以下为原代码

for (i in 1:nrow(file)){
reads=as.numeric(file[i,])
ref=reads[1:(ncol(file)/2)]
var=reads[(ncol(file)/2+1):ncol(file)]
tid=rep(1:(length(sel)/2),each=2)-1
sl=which(!is.na(ref))
if ((length(sl)/2)>=1){
if ((length(sl)/2)>=2){
#df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x1=sts[sl],x2=tid[sl])
df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x2=tid[sl])
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1 + f(x2, model = "iid")
#formula = ALT_COUNT ~ 1 + f(TISSUE_ID, model = "iid")
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}else{
#df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl],x1=sts[sl])
df=data.frame(y=var[sl],Ntrials=ref[sl]+var[sl])
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}
}
}